Skip to content

Speed up LinearQubitOperator matvec - #1417

Merged
mhucka merged 3 commits into
quantumlib:mainfrom
blazingphoenix7:vectorize-linear-qubit-operator-matvec
Jul 17, 2026
Merged

mhucka merged 3 commits into
quantumlib:mainfrom
blazingphoenix7:vectorize-linear-qubit-operator-matvec

Conversation

@blazingphoenix7

Copy link
Copy Markdown
Contributor

Second of the operator hot-path findings from #1407; the first is #1415.

LinearQubitOperator._matvec applies a QubitOperator to a state vector, and it
is the inner loop behind generate_linear_qubit_operator and the eigensolvers
built on top of it. For each term it took the input vector, split it into 2**k
pieces to reach the first non-identity qubit, split each of those in half, applied
the 2x2 Pauli to the pairs, and concatenated everything back, once per Pauli in
the term. That is a lot of Python-level list building and array copying for what
is, per term, a single reindexing of the amplitudes.

A Pauli string acts on the state vector as a signed permutation, so it can be
written directly on the flat index. Number the amplitudes 0 .. 2**n - 1 and put
qubit q at bit n - 1 - q (the ordering the split version already used, and the
one qubit_operator_sparse uses). Then for one term:

  • the X and Y qubits flip their bit, so amplitude c moves to c ^ x_mask;
  • the Y and Z qubits give a factor of -1 whenever the matching bit of c is set,
    that is (-1) raised to the parity of popcount(c & z_mask);
  • each Y gives an extra factor of 1j.

So the term becomes out[idx ^ x_mask] += coeff * 1j**(y_count % 4) * sign * x,
where sign is a vector of +-1 from the bit parity and idx ^ x_mask is a
permutation of the indices, which is what makes the fancy-indexed += safe. The
masks are cheap integer folds over the term and the per-amplitude work is a few
vectorized numpy operations with no intermediate lists.

This is an implementation change only. Public behavior is unchanged: complex
output, inputs of shape (dim,) or (dim, 1), the identity term, operators
declared on more qubits than they act on, and the empty operator all behave as
before. Only the serial _matvec is touched; ParallelLinearQubitOperator is
left alone.

The parity of popcount uses a small xor-fold helper on uint64 rather than
numpy.bitwise_count, because the project still supports numpy 1.26 (the
max_compat CI job pins numpy==1.26.4) and bitwise_count needs numpy 2.0.

Correctness

linear_qubit_operator_test.py passes unchanged. On top of that I compared the
new matvec against qubit_operator_sparse(op) @ x over random operators for
n = 2..12 (mixed X/Y/Z, term counts including the empty and identity-only
operators), real and complex x, shapes (dim,) and (dim, 1), and n_qubits
oversized by 2, plus a scipy.sparse.linalg.eigsh run checked against dense
eigenvalues. Worst absolute differences:

random real x, shape (dim,)          5.4e-15
random complex x, shape (dim,)       9.1e-15
random complex x, shape (dim, 1)     9.1e-15
oversized n_qubits by 2              2.4e-15
identity-only and empty operator     0
eigsh vs dense (3 smallest)          8.9e-15

black and mypy are clean on the file.

Motivation

This is what #499 asked for. That issue called for a better _matvec, measured
that one was worth 80-100x at system sizes of 16 qubits and up, and noted that
ParallelLinearQubitOperator could then be deleted as obsolete once the serial
path was fast. The sketch there routed each Pauli through Cirq's
apply_unitary; no port ever landed and the issue was closed for inactivity in
2025. This PR gets the same win with plain numpy on the flat index, so it adds
no dependency, and the 120x at 16 qubits below bears out the issue's estimate.

The parallel path is also where the deadlock possibility in #1405 lives, and
with the serial matvec this fast there is rarely a reason to reach for it.
Removing it, as #499 suggested, would be a separate change; this PR fixes
neither issue.

Benchmarks

Random QubitOperators with mixed X/Y/Z terms applied to a complex state vector:

n_qubits   terms    before     after    speedup
   16       100     20.5 s     0.17 s      120x
   18       100     87.0 s     1.5 s        60x
   20        50      157 s     3.0 s        52x

Measured on my laptop (Intel Core Ultra 5 225, Windows, Python 3.12, numpy 2.5).
I ran the base/patch comparison twice in alternation with fixed seeds, taking
medians over three runs at n = 16 and single runs at n = 18 and 20 since the old
code takes minutes there; the per-round speedups agreed within a few percent, so
the table pools both rounds.

_matvec applied each Pauli term by recursively splitting the amplitude
vector into sublists and reassembling it. Apply each term instead as a
signed permutation of the amplitudes: the X and Y qubits give an index
xor mask, the Y and Z qubits give the indices whose bit parity sets the
sign, and each Y adds a factor of 1j. This is a pure implementation
change with identical results.

@gemini-code-assist gemini-code-assist Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Code Review

This pull request refactors the _matvec method in LinearQubitOperator to use bitwise operations and masks, introducing a new _bit_parity helper function to apply operators as signed permutations instead of recursively splitting vectors. The review feedback points out that using 1j ** (y_count % 4) introduces numerical inaccuracies due to floating-point approximations and suggests using a lookup table [1, 1j, -1, -1j][y_count % 4] to keep the real and imaginary parts exact.

Comment thread src/openfermion/linalg/linear_qubit_operator.py Outdated
The phase contributed by the Y factors is always one of 1, 1j, -1, -1j.
Index a small table by y_count % 4 instead of evaluating 1j ** (y_count % 4).
Same values, and it stays exact if the exponent is ever a numpy integer
rather than a Python int.

@mhucka mhucka left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for this nice work! I had just a few suggestions.

Comment thread src/openfermion/linalg/linear_qubit_operator.py Outdated
Comment thread src/openfermion/linalg/linear_qubit_operator.py Outdated
Comment thread src/openfermion/linalg/linear_qubit_operator.py Outdated
@mhucka mhucka added the area/performance Involves code performance label Jul 16, 2026
Use a tuple rather than a list for the (1, 1j, -1, -1j) phase constant so
it is stored once instead of rebuilt per call, keep the sign array as int8
since its values are only +1 and -1, and allocate the index array lazily so
an operator with only the identity term skips it.

@mhucka mhucka left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the updates. Looks good.

@mhucka
mhucka added this pull request to the merge queue Jul 17, 2026
Merged via the queue into quantumlib:main with commit e8729f7 Jul 17, 2026
23 checks passed
mhucka pushed a commit to mhucka/OpenFermion that referenced this pull request Jul 27, 2026
…ransform (quantumlib#1427)

Third of the operator hot-path findings from quantumlib#1407; the earlier two are
quantumlib#1415 and quantumlib#1417.

`jordan_wigner` on an `InteractionOperator` maps the one- and two-body
integrals
straight to Pauli terms instead of going through a `FermionOperator`
first, and
quantumlib#587 points at exactly this routine as the fast path to build on ("there
is a
special routine that maps the InteractionOperators directly to qubits
... and
that is faster"). quantumlib#121 is the other side of it: it flagged this same
routine as
"surprisingly slow" back in 2017 and got closed without a fix.

The transform walks every one-body pair and every pair of pairs, so the
two-body
part is O(N^4). For each Pauli string it built a fresh single-term
`QubitOperator` and summed that into the running total. Each of those
constructions parses the term, validates every `(index, action)` factor,
runs the
Pauli `_simplify`, and in the three-unique-index case multiplies by a
`Z`
operator, which deepcopies. None of that is needed to add a term to a
sum: a
`QubitOperator`'s `terms` is just a dict from a sorted tuple of `(index,
action)`
factors to a coefficient, and the strings this routine emits are already
sorted
with their factors already combined, so there is nothing for the
constructor to
normalize.

This accumulates the `pauli_string -> coefficient` pairs into a plain
dict with a
small `_add_term` helper that does what `__iadd__` / `__isub__` do to a
single
term (add onto the running value, drop it once it falls under
`EQ_TOLERANCE`), and
builds one `QubitOperator` at the end by assigning the finished dict to
`.terms`.
Accumulating into a dict instead of building and summing single-term
operators
avoids that per-term construction and validation, and it measures about
2x
faster, with identical output.

### Public functions

`jordan_wigner_one_body` and `jordan_wigner_two_body` are public and
re-exported,
and `bksf_test` and `jordan_wigner_test` call them directly, so they
keep their
signatures and their output. I pulled the term-building logic into two
private
helpers (`_one_body_terms`, `_two_body_terms`) that return the dict, and
the
public functions and the interaction-op driver both go through them, so
there is
one copy of the logic and the public results are unchanged. The `Z`
times
hopping-term product in the three-unique case is a single-qubit
relabeling (the
shared qubit carries at most a parity `Z`, and `Z` squared is the
identity), so it
is a direct toggle of that index in the string now rather than an
operator
multiply and deepcopy.

### Correctness

Before editing the source I captured the current `jordan_wigner` output
on 104
random real `InteractionOperator`s and saved it. The mix is physical
(symmetric, molecular-style) operators, dense generic real tensors,
sparse and
very sparse tensors, one-body-only, constant-only, and the exact-zero
operator,
at N = 1, 2, 3, 4, 6, 8, 10, 12. The new code reproduces that captured
output
exactly on every case: identical term-key sets, max absolute coefficient
difference 0, and `(new - old).induced_norm(1)` equal to 0. It also
matches an
independent dict-based reference I wrote separately, to the same 0.
Seeds are
fixed, so the check is deterministic.

`jordan_wigner_test.py` and `bksf_test.py` pass (the two public helpers
are
checked there against the `FermionOperator` path over a full index grid
with
complex coefficients), and `transforms/` passes as a whole. black and
mypy are
clean on the file. Every changed line is hit by the existing suite
except the
pre-existing `n_qubits < count_qubits(iop)` guard, which the old code
did not
cover either and which I left untouched.

The routine is still documented as real-only and this does not change
that. It
computes the same coefficients from the same integrals; only the
container they
land in changed.

### Benchmarks

Random real InteractionOperators, cumulative `jordan_wigner` time.
Medians over 5
repeats within each of 6 alternated upstream/branch rounds, fixed seeds,
on an
Intel Core Ultra 5 225 laptop (Windows, Python 3.12, numpy 2.5).

    N     upstream    ours       speedup
    8     16 ms       7.6 ms     2.1x
    12    101 ms      48 ms      2.1x
    16    359 ms      177 ms     2.0x
    20    963 ms      476 ms     2.0x
    24    2222 ms     1093 ms    2.0x
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

area/performance Involves code performance

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants